{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# dislocation_monopole calculation style\n",
    "\n",
    "**Lucas M. Hale**, [lucas.hale@nist.gov](mailto:lucas.hale@nist.gov?Subject=ipr-demo), *Materials Science and Engineering Division, NIST*.\n",
    "\n",
    "## Introduction\n",
    "\n",
    "The dislocation_monopole calculation style calculation inserts a dislocation monopole into an atomic system using the anisotropic elasticity solution for a perfectly straight dislocation. The system is divided into two regions: a boundary region at the system's edges perpendicular to the dislocation line, and an active region around the dislocation. After inserting the dislocation, the atoms within the active region of the system are relaxed while the positions of the atoms in the boundary region are held fixed at the elasticity solution. The relaxed dislocation system and corresponding dislocation-free base systems are retained in the calculation's archived record. Various properties associated with the dislocation's elasticity solution are recorded in the calculation's results record.\n",
    "\n",
    "### Version notes\n",
    "\n",
    "- 2018-09-25: Notebook added\n",
    "- 2019-07-30: Method and Notebook updated for iprPy version 0.9.\n",
    "- 2020-09-22: Notebook updated to reflect that calculation method has changed to now use atomman.defect.Dislocation.  \n",
    "\n",
    "### Additional dependencies\n",
    "\n",
    "### Disclaimers\n",
    "\n",
    "- [NIST disclaimers](http://www.nist.gov/public_affairs/disclaimer.cfm)\n",
    "- This calculation method holds the boundary atoms fixed during the relaxation process which results in a mismatch stress at the border between the active and fixed regions that interacts with the dislocation.  Increasing the distance between the dislocation and the boundary region, i.e. increasing the system size, will reduce the influence of the mismatch stresses.\n",
    "- The boundary atoms are fixed at the elasticity solution, which assumes the dislocation to be compact (not spread out, dissociated into partials) and to remain at the center of the system.  An alternate simulation design or boundary region method should be used if you want accurate simulations of dislocations with wide cores and/or in which the dislocation moves long distances.\n",
    "- The calculation allows for the system to be relaxed either using only static energy/force minimizations or with molecular dynamic steps followed by a minimization.  Only performing a static relaxation is considerably faster than performing a dynamic relaxation, but the dynamic relaxation is more likely to result in a better final dislocation structure.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method and Theory\n",
    "\n",
    "### Stroh theory\n",
    "\n",
    "A detailed description of the Stroh method to solve the Eshelby solution for an anisotropic straight dislocation can be found in the [atomman documentation](https://www.ctcms.nist.gov/potentials/atomman/).\n",
    "\n",
    "### Orientation\n",
    "\n",
    "One of the most important considerations in defining an atomistic system containing a dislocation monopole system is the system's orientation.  In particular, care is needed to properly align the system's Cartesian axes, $x, y, z$, the system's box vectors, $a, b, c$, and the Stroh solution's axes, $u, m, n$.\n",
    "\n",
    "- In order for the dislocation to remain perfectly straight in the atomic system, the dislocation line, $u$, must correspond to one of the system's box vectors.  The resulting dislocation monopole system will be periodic along the box direction corresponding to $u$, and non-periodic in the other two box directions.\n",
    "\n",
    "- The Stroh solution is invariant along the dislocation line direction, $u$, thereby the solution is 2 dimensional. $m$ and $n$ are the unit vectors for the 2D axis system used by the Stroh solution. As such, $u, m$ and $n$ are all normal to each other. The plane defined by the $um$ vectors is the dislocation's slip plane, i.e. $n$ is normal to the slip plane.\n",
    "\n",
    "- For LAMMPS simulations, the system's box vectors are limited such that $a$ is parallel to the $x$-axis, and $b$ is in the $xy$-plane (i.e. has no $z$ component). Based on this and the previous two points, the most convenient, and therefore the default, orientation for a generic dislocation is to\n",
    "\n",
    "  - Make $u$ and $a$ parallel, which also places $u$ parallel to the $x$-axis.\n",
    "\n",
    "  - Select $b$ such that it is within the slip plane. As $u$ and $a$ must also be in the slip plane, the plane itself is defined by $a \\times b$.\n",
    "  \n",
    "  - Given that neither $a$ nor $b$ have $z$ components, the normal to the slip plane will be perpendicular to $z$.  From this, it naturally follows that $m$ can be taken as parallel to the $y$-axis, and $n$ parallel to the $z$-axis.\n",
    "\n",
    "### Calculation methodology\n",
    "\n",
    "1. An initial system is generated based on the loaded system and *uvw*, *atomshift*, and *sizemults* input parameters.  This initial system is saved as base.dump.\n",
    "\n",
    "2. The atomman.defect.Stroh class is used to obtain the dislocation solution based on the defined $m$ and $n$ vectors, $C_{ij}$, and the dislocation's Burgers vector, $b_i$.\n",
    "\n",
    "3. The dislocation is inserted into the system by displacing all atoms according to the Stroh solution's displacements.  The dislocation line is placed parallel to the specified *dislocation_lineboxvector* and includes the Cartesian point (0, 0, 0).\n",
    "\n",
    "4. The box dimension parallel to the dislocation line is kept periodic, and the other two box directions are made non-periodic. A boundary region is defined that is at least *bwidth* thick at the edges of the two non-periodic directions, such that the shape of the active region corresponds to the *bshape* parameter. Atoms in the boundary region are identified by altering their integer atomic types.\n",
    "\n",
    "5. The dislocation is relaxed by performing a LAMMPS simulation.  The atoms in the active region are allowed to relax either with nvt integration followed by an energy/force minimization, or with just an energy/force minimization.  The atoms in the boundary region are held fixed at the elastic solution.  The resulting relaxed system is saved as disl.dump.\n",
    "\n",
    "6. Parameters associated with the Stroh solution are saved to the results record.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Demonstration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.1. Library imports"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import libraries needed by the calculation. The external libraries used are:\n",
    "\n",
    "- [numpy](http://www.numpy.org/)\n",
    "\n",
    "- [atomman](https://github.com/usnistgov/atomman)\n",
    "\n",
    "- [iprPy](https://github.com/usnistgov/iprPy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Notebook last executed on 2020-09-22 using iprPy version 0.10.2\n"
     ]
    }
   ],
   "source": [
    "# Standard library imports\n",
    "from pathlib import Path\n",
    "import os\n",
    "import random\n",
    "import datetime\n",
    "\n",
    "# http://www.numpy.org/\n",
    "import numpy as np \n",
    "\n",
    "# https://github.com/usnistgov/atomman \n",
    "import atomman as am\n",
    "import atomman.lammps as lmp\n",
    "import atomman.unitconvert as uc\n",
    "\n",
    "# https://github.com/usnistgov/iprPy\n",
    "import iprPy\n",
    "\n",
    "print('Notebook last executed on', datetime.date.today(), 'using iprPy version', iprPy.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 3, 6])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = np.array([1,3,-6])\n",
    "np.abs(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.2. Default calculation setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify calculation style\n",
    "calc_style = 'dislocation_monopole'\n",
    "\n",
    "# If workingdir is already set, then do nothing (already in correct folder)\n",
    "try:\n",
    "    workingdir = workingdir\n",
    "\n",
    "# Change to workingdir if not already there\n",
    "except:\n",
    "    workingdir = Path('calculationfiles', calc_style)\n",
    "    if not workingdir.is_dir():\n",
    "        workingdir.mkdir(parents=True)\n",
    "    os.chdir(workingdir)\n",
    "    \n",
    "# Initialize connection to library\n",
    "with open('C:/Users/lmh1/Documents/potentials_nist_gov/password.txt') as f:\n",
    "    user, pswd = f.read().strip().split()\n",
    "library = iprPy.Library(load=['lammps_potentials'], username=user, password=pswd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Assign values for the calculation's run parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1. Specify system-specific paths\n",
    "\n",
    "- __lammps_command__ (required) is the LAMMPS command to use.\n",
    "\n",
    "- __mpi_command__ MPI command for running LAMMPS in parallel. A value of None will run simulations serially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "lammps_command = 'lmp_serial'\n",
    "mpi_command = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2. Load interatomic potential\n",
    "\n",
    "- __potential_name__ gives the name of the potential_LAMMPS reference record in the iprPy library to use for the calculation.  \n",
    "\n",
    "- __potential__ is an atomman.lammps.Potential object (required)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "potential_name = '1999--Mishin-Y--Ni--LAMMPS--ipr1'\n",
    "\n",
    "# Retrieve potential and parameter file(s)\n",
    "potential = library.get_lammps_potential(id=potential_name, getfiles=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.3. Load initial unit cell system\n",
    "\n",
    "- __ucell__ is an atomman.System representing a fundamental unit cell of the system (required).  Here, this is loaded from the database for the prototype."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [ 3.520,  0.000,  0.000]\n",
      "bvect =  [ 0.000,  3.520,  0.000]\n",
      "cvect =  [ 0.000,  0.000,  3.520]\n",
      "origin = [ 0.000,  0.000,  0.000]\n",
      "natoms = 4\n",
      "natypes = 1\n",
      "symbols = ('Ni',)\n",
      "pbc = [ True  True  True]\n",
      "per-atom properties = ['atype', 'pos']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |   0.000 |   0.000 |   0.000\n",
      "      1 |       1 |   0.000 |   1.760 |   1.760\n",
      "      2 |       1 |   1.760 |   0.000 |   1.760\n",
      "      3 |       1 |   1.760 |   1.760 |   0.000\n"
     ]
    }
   ],
   "source": [
    "# Create ucell by loading prototype record\n",
    "ucell = am.load('crystal', potential=potential, family='A1--Cu--fcc', database=library)\n",
    "\n",
    "print(ucell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.4 Specify material elastic constants\n",
    "\n",
    "Simple input parameters:\n",
    "\n",
    "- __C_dict__ is a dictionary containing the unique elastic constants for the potential and crystal structure defined above. \n",
    "\n",
    "Derived parameters\n",
    "\n",
    "- __C__ is an atomman.ElasticConstants object built from C_dict."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "C_dict = {}\n",
    "C_dict['C11'] = uc.set_in_units(247.86, 'GPa')\n",
    "C_dict['C12'] = uc.set_in_units(147.83, 'GPa')\n",
    "C_dict['C44'] = uc.set_in_units(124.84, 'GPa')\n",
    "\n",
    "# -------------- Derived parameters -------------- #\n",
    "# Build ElasticConstants object from C_dict terms\n",
    "C = am.ElasticConstants(**C_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.5 Specify the defect parameters\n",
    "\n",
    "- __burgers__ is the crystallographic Miller Burgers vector for the dislocation. \n",
    "\n",
    "- __ξ_uvw__ is the Miller \\[uvw\\] line vector direction for the dislocation.  The angle between burgers and ξ_uvw determines the dislocation's character\n",
    "\n",
    "- __slip_hkl__ is the Miller (hkl) slip plane for the dislocation.\n",
    "    \n",
    "- __m__ is the Cartesian vector of the final system that the dislocation solution's m vector (in-plane, perpendicular to ξ) should align with.  Limited to being parallel to one of the three Cartesian axes.  \n",
    "\n",
    "- __n__ is the Cartesian vector of the final system that the dislocation solution's n vector (slip plane normal) should align with.  Limited to being parallel to one of the three Cartesian axes. \n",
    "\n",
    "- __shift__ is a rigid body shift to apply to the atoms in the system. This controls how the atomic positions align with the ideal position of the dislocation core, which is at coordinates (0,0) for the two Cartesian axes aligned with m and n.\n",
    "\n",
    "- __shiftscale__ allows for shift to be defined relative to the cell created by rotating ucell to coincide with the dislocation solution orientation.  This is useful as it allows for shift values to be defined relative to the defect type and crystal prototype rather than on a per-crystal basis. \n",
    "\n",
    "- __shiftindex__ alternate to specifying shift values, the shiftindex allows for one of the identified suggested shift values to be used that will position the slip plane halfway between two planes of atoms.  Note that shiftindex values only shift atoms in the slip plane normal direction and may not be the ideal positions for some dislocation cores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fcc a/2 <110>{111} dislocations\n",
    "burgers = 0.5 * np.array([ 1., -1., 0.])\n",
    "slip_hkl = [1, 1, 1]\n",
    "\n",
    "# Line direction determines dislocation character\n",
    "ξ_uvw = [ 1, 1,-2] # 90 degree edge\n",
    "#ξ_uvw = [ 1, 0,-1] # 60 degree mixed\n",
    "#ξ_uvw = [ 1,-2, 1] # 30 degree mixed\n",
    "#ξ_uvw = [ 1,-1, 0] # 0 degree screw\n",
    "\n",
    "# Best choice for m + n as it works for non-cubic systems\n",
    "m = [0,1,0]\n",
    "n = [0,0,1]\n",
    "\n",
    "# Specify shift or shiftindex\n",
    "shift = None\n",
    "shiftscale = True\n",
    "shiftindex = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.6 Specify calculation-specific run parameters\n",
    "\n",
    "Simple input parameters:\n",
    "\n",
    "- __boundarywidth__ sets the minimum width of the fixed-atom boundary region.\n",
    "\n",
    "- __boundaryscale__ flag indicating if boundarywidth is absolute (False) or scaled relative to ucell's a lattice parameter (True).\n",
    "\n",
    "- __boundaryshape__ specifies what shape to make the fixed-atom boundary region.\n",
    "\n",
    "    - __cylinder__ will create a cylindrical active region.\n",
    "    \n",
    "    - __box__ will create a rectangular active region.\n",
    "\n",
    "- __annealtemperature__ is the temperature at which to relax (anneal) the dislocation system. If annealtemperature is 0.0, then only a static relaxation will be performed. Default value is 0.0.\n",
    "\n",
    "- __annealsteps__ is the number of nvt iteration steps to perform at the given temperature.  Default value is 0 if annealtemperature is zero, 10000 otherwise.\n",
    "\n",
    "- __randomseed__ allows for the random seed used in generating initial atomic velocities for a dynamic relaxation to be specified. This is an integer between 1 and 900000000. Default value is None, which will randomly pick a number in that range.\n",
    "\n",
    "- __energytolerance__ is the energy tolerance to use during the minimizations. This is unitless.\n",
    "\n",
    "- __forcetolerance__ is the force tolerance to use during the minimizations. This is in energy/length units.\n",
    "\n",
    "- __maxiterations__ is the maximum number of minimization iterations to use.\n",
    "\n",
    "- __maxevaluations__ is the maximum number of minimization evaluations to use.\n",
    "\n",
    "- __maxatommotion__ is the largest distance that an atom is allowed to move during a minimization iteration. This is in length units. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify the boundary shape to use\n",
    "boundarywidth = 3\n",
    "boundaryscale = True\n",
    "boundaryshape = 'cylinder'\n",
    "\n",
    "# Specify MD anneal parameters\n",
    "annealtemperature = 50.0\n",
    "annealsteps = 10000\n",
    "randomseed = None\n",
    "\n",
    "# Specify minimization parameters\n",
    "energytolerance = 0.0\n",
    "forcetolerance = uc.set_in_units(1e-6, 'eV/angstrom')\n",
    "maxiterations = 10000\n",
    "maxevaluations = 100000\n",
    "maxatommotion = uc.set_in_units(0.01, 'angstrom')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.7. Modify system\n",
    "\n",
    "- __sizemults__ list of three sets of two integers specifying how many times the ucell vectors of $a$, $b$ and $c$ are replicated in positive and negative directions when creating system.\n",
    "\n",
    "- __amin, bmin, cmin__ specifies minimum widths in the three directions. The corresponding sizemults will be increased to ensure these widths if needed.\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "sizemults = [1, 40, 40]\n",
    "amin = uc.set_in_units(0.0, 'angstrom')\n",
    "bmin = uc.set_in_units(0.0, 'angstrom')\n",
    "cmin = uc.set_in_units(0.0, 'angstrom')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Define calculation function(s) and generate template LAMMPS script(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1. disl_relax.template"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('disl_relax.template', 'w') as f:\n",
    "    f.write(\"\"\"#LAMMPS input script for relaxing a dislocation\n",
    "\n",
    "<atomman_system_pair_info>\n",
    "\n",
    "group move type <group_move>\n",
    "group hold subtract all move\n",
    "\n",
    "compute peatom all pe/atom\n",
    "\n",
    "dump first all custom <maxeval> *.dump id type x y z c_peatom\n",
    "dump_modify first format <dump_modify_format>\n",
    "thermo_style custom step pe\n",
    "\n",
    "fix nomove hold setforce 0.0 0.0 0.0\n",
    "\n",
    "<anneal_info>\n",
    "\n",
    "min_modify dmax <dmax>\n",
    "\n",
    "minimize <etol> <ftol> <maxiter> <maxeval>\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2. dislocationmonopole()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dislocationmonopole(lammps_command, ucell, potential, C, burgers, ξ_uvw,\n",
    "                        slip_hkl, mpi_command=None, m=[0,1,0], n=[0,0,1],\n",
    "                        sizemults=None, amin=None, bmin=None, cmin=None,\n",
    "                        shift=None, shiftscale=False, shiftindex=None, tol=1e-8,\n",
    "                        etol=0.0, ftol=0.0, maxiter=10000, maxeval=100000,\n",
    "                        dmax=uc.set_in_units(0.01, 'angstrom'),\n",
    "                        annealtemp=0.0, annealsteps=None, randomseed=None,\n",
    "                        boundaryshape='cylinder', boundarywidth=0.0, boundaryscale=False):\n",
    "    \"\"\"\n",
    "    Creates and relaxes a dislocation monopole system.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    lammps_command :str\n",
    "        Command for running LAMMPS.\n",
    "    ucell : atomman.System\n",
    "        The unit cell to use as the seed for generating the dislocation\n",
    "        monopole system.\n",
    "    potential : atomman.lammps.Potential\n",
    "        The LAMMPS implemented potential to use.\n",
    "    C : atomman.ElasticConstants\n",
    "        The elastic constants associated with the bulk crystal structure\n",
    "        for ucell.\n",
    "    burgers : array-like object\n",
    "        The dislocation's Burgers vector given as a Miller or\n",
    "        Miller-Bravais vector relative to ucell.\n",
    "    ξ_uvw : array-like object\n",
    "        The dislocation's line direction given as a Miller or\n",
    "        Miller-Bravais vector relative to ucell.\n",
    "    slip_hkl : array-like object\n",
    "        The dislocation's slip plane given as a Miller or Miller-Bravais\n",
    "        plane relative to ucell.\n",
    "    mpi_command : str or None, optional\n",
    "        The MPI command for running LAMMPS in parallel.  If not given, LAMMPS\n",
    "        will run serially.\n",
    "    m : array-like object, optional\n",
    "        The m unit vector for the dislocation solution.  m, n, and ξ\n",
    "        (dislocation line) should be right-hand orthogonal.  Default value\n",
    "        is [0,1,0] (y-axis).\n",
    "    n : array-like object, optional\n",
    "        The n unit vector for the dislocation solution.  m, n, and ξ\n",
    "        (dislocation line) should be right-hand orthogonal.  Default value\n",
    "        is [0,0,1] (z-axis). n is normal to the dislocation slip plane.\n",
    "    sizemults : tuple, optional\n",
    "        The size multipliers to use when generating the system.  Values are\n",
    "        limited to being positive integers.  The multipliers for the two\n",
    "        non-periodic directions must be even.  If not given, the default\n",
    "        multipliers will be 2 for the non-periodic directions and 1 for the\n",
    "        periodic direction.\n",
    "    amin : float, optional\n",
    "        A minimum thickness to use for the a box vector direction of the\n",
    "        final system.  Default value is 0.0.  For the non-periodic\n",
    "        directions, the resulting vector multiplier will be even.  If both\n",
    "        amin and sizemults is given, then the larger multiplier for the two\n",
    "        will be used.\n",
    "    bmin : float, optional\n",
    "        A minimum thickness to use for the b box vector direction of the\n",
    "        final system.  Default value is 0.0.  For the non-periodic\n",
    "        directions, the resulting vector multiplier will be even.  If both\n",
    "        bmin and sizemults is given, then the larger multiplier for the two\n",
    "        will be used.\n",
    "    cmin : float, optional\n",
    "        A minimum thickness to use for the c box vector direction of the\n",
    "        final system.  Default value is 0.0.  For the non-periodic\n",
    "        directions, the resulting vector multiplier will be even.  If both\n",
    "        cmin and sizemults is given, then the larger multiplier for the two\n",
    "        will be used.\n",
    "    shift : float, optional\n",
    "        A rigid body shift to apply to the rotated cell prior to inserting\n",
    "        the dislocation.  Should be selected such that the ideal slip plane\n",
    "        does not correspond to any atomic planes.  Is taken as absolute if\n",
    "        shiftscale is False, or relative to the rotated cell's box vectors\n",
    "        if shiftscale is True.  Cannot be given with shiftindex.  If\n",
    "        neither shift nor shiftindex is given then shiftindex = 0 is used.\n",
    "    shiftindex : float, optional\n",
    "        The index of the identified optimum shifts based on the rotated\n",
    "        cell to use.  Different values allow for the selection of different\n",
    "        atomic planes neighboring the slip plane.  Note that shiftindex\n",
    "        values only apply shifts normal to the slip plane; best shifts for\n",
    "        non-planar dislocations (like bcc screw) may also need a shift in\n",
    "        the slip plane.  Cannot be given with shiftindex.  If neither shift\n",
    "        nor shiftindex is given then shiftindex = 0 is used.\n",
    "    shiftscale : bool, optional\n",
    "        If False (default), a given shift value will be taken as absolute\n",
    "        Cartesian.  If True, a given shift will be taken relative to the\n",
    "        rotated cell's box vectors.\n",
    "    tol : float\n",
    "        A cutoff tolerance used with obtaining the dislocation solution.\n",
    "        Only needs to be changed if there are issues with obtaining a\n",
    "        solution.\n",
    "    etol : float, optional\n",
    "        The energy tolerance for the structure minimization. This value is\n",
    "        unitless. Default is 0.0.\n",
    "    ftol : float, optional\n",
    "        The force tolerance for the structure minimization. This value is in\n",
    "        units of force. Default is 0.0.\n",
    "    maxiter : int, optional\n",
    "        The maximum number of minimization iterations to use. Default is \n",
    "        10000.\n",
    "    maxeval : int, optional\n",
    "        The maximum number of minimization evaluations to use. Default is \n",
    "        100000.\n",
    "    dmax : float, optional\n",
    "        The maximum distance in length units that any atom is allowed to relax\n",
    "        in any direction during a single minimization iteration. Default is\n",
    "        0.01 Angstroms.\n",
    "    annealtemp : float, optional\n",
    "        The temperature to perform a dynamic relaxation at. Default is 0.0,\n",
    "        which will skip the dynamic relaxation.\n",
    "    annealsteps : int, optional\n",
    "        The number of time steps to run the dynamic relaxation for.  Default\n",
    "        is None, which will run for 10000 steps if annealtemp is not 0.0.  \n",
    "    randomseed : int or None, optional\n",
    "        Random number seed used by LAMMPS in creating velocities and with\n",
    "        the Langevin thermostat.  Default is None which will select a\n",
    "        random int between 1 and 900000000.\n",
    "    boundaryshape : str, optional\n",
    "        Indicates the shape of the boundary region to use.  Options are\n",
    "        'cylinder' (default) and 'box'.  For 'cylinder', the non-boundary\n",
    "        region is defined by a cylinder with axis along the dislocation\n",
    "        line and a radius that ensures the boundary is at least\n",
    "        boundarywidth thick.  For 'box', the boundary region will be\n",
    "        exactly boundarywidth thick all around.      \n",
    "    boundarywidth : float, optional\n",
    "        The width of the boundary region to apply.  Default value is 0.0,\n",
    "        i.e. no boundary region.  All atoms in the boundary region will\n",
    "        have their atype values changed.\n",
    "    boundaryscale : bool, optional\n",
    "        If False (Default), the boundarywidth will be taken as absolute.\n",
    "        If True, the boundarywidth will be taken relative to the magnitude\n",
    "        of the unit cell's a box vector.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    dict\n",
    "        Dictionary of results consisting of keys:\n",
    "        \n",
    "        - **'dumpfile_base'** (*str*) - The filename of the LAMMPS dump file\n",
    "          for the relaxed base system.\n",
    "        - **'symbols_base'** (*list of str*) - The list of element-model\n",
    "          symbols for the Potential that correspond to the base system's\n",
    "          atypes.\n",
    "        - **'dumpfile_disl'** (*str*) - The filename of the LAMMPS dump file\n",
    "          for the relaxed dislocation monopole system.\n",
    "        - **'symbols_disl'** (*list of str*) - The list of element-model\n",
    "          symbols for the Potential that correspond to the dislocation\n",
    "          monopole system's atypes.\n",
    "        - **'dislocation'** (*atomman.defect.Dislocation*) - The Dislocation\n",
    "          object used to generate the monopole system.\n",
    "        - **'E_total_disl'** (*float*) - The total potential energy of the\n",
    "          dislocation monopole system.\n",
    "    \"\"\"\n",
    "    # Construct dislocation configuration generator\n",
    "    dislocation = am.defect.Dislocation(ucell, C, burgers, ξ_uvw, slip_hkl,\n",
    "                                        m=m, n=n, shift=shift, shiftindex=shiftindex,\n",
    "                                        shiftscale=shiftscale, tol=tol)\n",
    "    \n",
    "    # Generate the base and dislocation systems\n",
    "    base_system, disl_system = dislocation.monopole(sizemults=sizemults,\n",
    "                                                    amin=amin, bmin=bmin, cmin=cmin,\n",
    "                                                    shift=shift,\n",
    "                                                    shiftindex=shiftindex,\n",
    "                                                    shiftscale=shiftscale,\n",
    "                                                    boundaryshape=boundaryshape,\n",
    "                                                    boundarywidth=boundarywidth,\n",
    "                                                    boundaryscale=boundaryscale,\n",
    "                                                    return_base_system=True)\n",
    "    \n",
    "    # Initialize results dict\n",
    "    results_dict = {}\n",
    "    \n",
    "    # Save initial perfect system\n",
    "    base_system.dump('atom_dump', f='base.dump')\n",
    "    results_dict['dumpfile_base'] = 'base.dump'\n",
    "    results_dict['symbols_base'] = base_system.symbols\n",
    "    \n",
    "    # Save dislocation generator\n",
    "    results_dict['dislocation'] = dislocation\n",
    "\n",
    "    # Relax system\n",
    "    relaxed = disl_relax(lammps_command, disl_system, potential,\n",
    "                         mpi_command=mpi_command, annealtemp=annealtemp,\n",
    "                         annealsteps=annealsteps, randomseed=randomseed,\n",
    "                         etol=etol, ftol=ftol, maxiter=maxiter, \n",
    "                         maxeval=maxeval, dmax=dmax)\n",
    "    \n",
    "    # Save relaxed dislocation system with original box vects\n",
    "    system_disl = am.load('atom_dump', relaxed['dumpfile'], symbols=disl_system.symbols)\n",
    "    \n",
    "    system_disl.box_set(vects=disl_system.box.vects, origin=disl_system.box.origin)\n",
    "    system_disl.dump('atom_dump', f='disl.dump')\n",
    "    results_dict['dumpfile_disl'] = 'disl.dump'\n",
    "    results_dict['symbols_disl'] = system_disl.symbols\n",
    "    \n",
    "    results_dict['E_total_disl'] = relaxed['E_total']\n",
    "    \n",
    "    # Cleanup files\n",
    "    Path('0.dump').unlink()\n",
    "    Path(relaxed['dumpfile']).unlink()\n",
    "    for dumpjsonfile in Path('.').glob('*.dump.json'):\n",
    "        dumpjsonfile.unlink()\n",
    "    \n",
    "    return results_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.3. disl_relax()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def disl_relax(lammps_command, system, potential,\n",
    "               mpi_command=None, \n",
    "               annealtemp=0.0, annealsteps=None, randomseed=None,\n",
    "               etol=0.0, ftol=1e-6, maxiter=10000, maxeval=100000,\n",
    "               dmax=uc.set_in_units(0.01, 'angstrom')):\n",
    "    \"\"\"\n",
    "    Sets up and runs the disl_relax.in LAMMPS script for relaxing a\n",
    "    dislocation monopole system.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    lammps_command :str\n",
    "        Command for running LAMMPS.\n",
    "    system : atomman.System\n",
    "        The system to perform the calculation on.\n",
    "    potential : atomman.lammps.Potential\n",
    "        The LAMMPS implemented potential to use.\n",
    "    mpi_command : str, optional\n",
    "        The MPI command for running LAMMPS in parallel.  If not given, LAMMPS\n",
    "        will run serially.\n",
    "    annealtemp : float, optional\n",
    "        The temperature to perform a dynamic relaxation at. Default is 0.0,\n",
    "        which will skip the dynamic relaxation.\n",
    "    annealsteps : int, optional\n",
    "        The number of time steps to run the dynamic relaxation for.  Default\n",
    "        is None, which will run for 10000 steps if annealtemp is not 0.0.        \n",
    "    randomseed : int or None, optional\n",
    "        Random number seed used by LAMMPS in creating velocities and with\n",
    "        the Langevin thermostat.  Default is None which will select a\n",
    "        random int between 1 and 900000000.\n",
    "    etol : float, optional\n",
    "        The energy tolerance for the structure minimization. This value is\n",
    "        unitless. Default is 0.0.\n",
    "    ftol : float, optional\n",
    "        The force tolerance for the structure minimization. This value is in\n",
    "        units of force. Default is 0.0.\n",
    "    maxiter : int, optional\n",
    "        The maximum number of minimization iterations to use default is \n",
    "        10000.\n",
    "    maxeval : int, optional\n",
    "        The maximum number of minimization evaluations to use default is \n",
    "        100000.\n",
    "    dmax : float, optional\n",
    "        The maximum distance in length units that any atom is allowed to relax\n",
    "        in any direction during a single minimization iteration default is\n",
    "        0.01 Angstroms.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    dict\n",
    "        Dictionary of results consisting of keys:\n",
    "        \n",
    "        - **'logfile'** (*str*) - The name of the LAMMPS log file.\n",
    "        - **'dumpfile'** (*str*) - The name of the LAMMPS dump file\n",
    "          for the relaxed system.\n",
    "        - **'E_total'** (*float*) - The total potential energy for the\n",
    "          relaxed system.\n",
    "    \"\"\"\n",
    "    # Build filedict if function was called from iprPy\n",
    "    try:\n",
    "        assert __name__ == pkg_name\n",
    "        calc = iprPy.load_calculation(calculation_style)\n",
    "        filedict = calc.filedict\n",
    "    except:\n",
    "        filedict = {}\n",
    "\n",
    "    # Get lammps units\n",
    "    lammps_units = lmp.style.unit(potential.units)\n",
    "    \n",
    "    #Get lammps version date\n",
    "    lammps_date = lmp.checkversion(lammps_command)['date']\n",
    "    \n",
    "    # Define lammps variables\n",
    "    lammps_variables = {}\n",
    "    system_info = system.dump('atom_data', f='system.dat',\n",
    "                              potential=potential,\n",
    "                              return_pair_info=True)\n",
    "    lammps_variables['atomman_system_pair_info'] = system_info\n",
    "    lammps_variables['anneal_info'] = anneal_info(annealtemp, annealsteps, \n",
    "                                                  randomseed, potential.units)\n",
    "    lammps_variables['etol'] = etol\n",
    "    lammps_variables['ftol'] = uc.get_in_units(ftol, lammps_units['force'])\n",
    "    lammps_variables['maxiter'] = maxiter\n",
    "    lammps_variables['maxeval'] = maxeval\n",
    "    lammps_variables['dmax'] = dmax\n",
    "    lammps_variables['group_move'] = ' '.join(np.array(range(1, system.natypes // 2 + 1), dtype=str))\n",
    "    \n",
    "    # Set dump_modify format based on dump_modify_version\n",
    "    if lammps_date < datetime.date(2016, 8, 3):\n",
    "        lammps_variables['dump_modify_format'] = '\"%d %d %.13e %.13e %.13e %.13e\"'\n",
    "    else:\n",
    "        lammps_variables['dump_modify_format'] = 'float %.13e'\n",
    "    \n",
    "    # Write lammps input script\n",
    "    template_file = 'disl_relax.template'\n",
    "    lammps_script = 'disl_relax.in'\n",
    "    template = iprPy.tools.read_calc_file(template_file, filedict)\n",
    "    with open(lammps_script, 'w') as f:\n",
    "        f.write(iprPy.tools.filltemplate(template, lammps_variables,\n",
    "                                         '<', '>'))\n",
    "    \n",
    "    # Run LAMMPS\n",
    "    output = lmp.run(lammps_command, lammps_script, mpi_command)\n",
    "    thermo = output.simulations[-1]['thermo']\n",
    "    \n",
    "    # Extract output values\n",
    "    results = {}\n",
    "    results['logfile'] = 'log.lammps'\n",
    "    results['dumpfile'] = '%i.dump' % thermo.Step.values[-1] \n",
    "    results['E_total'] = uc.set_in_units(thermo.PotEng.values[-1],\n",
    "                                         lammps_units['energy'])\n",
    "    \n",
    "    return results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.4 anneal_info() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def anneal_info(temperature=0.0, runsteps=None, randomseed=None, units='metal'):\n",
    "    \"\"\"\n",
    "    Generates LAMMPS commands for thermo anneal.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    temperature : float, optional\n",
    "        The temperature to relax at (default is 0.0).\n",
    "    randomseed : int or None, optional\n",
    "        Random number seed used by LAMMPS in creating velocities and with\n",
    "        the Langevin thermostat.  (Default is None which will select a\n",
    "        random int between 1 and 900000000.)\n",
    "    units : str, optional\n",
    "        The LAMMPS units style to use (default is 'metal').\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    str\n",
    "        The generated LAMMPS input lines for performing a dynamic relax.\n",
    "        Will be '' if temperature==0.0.\n",
    "    \"\"\"    \n",
    "    # Return nothing if temperature is 0.0 (don't do thermo anneal)\n",
    "    if temperature == 0.0:\n",
    "        return ''\n",
    "    \n",
    "    # Generate velocity, fix nvt, and run LAMMPS command lines\n",
    "    else:\n",
    "        if randomseed is None:\n",
    "            randomseed = random.randint(1, 900000000)\n",
    "        if runsteps is None:\n",
    "            runsteps = 10000\n",
    "        \n",
    "        start_temp = 2 * temperature\n",
    "        tdamp = 100 * lmp.style.timestep(units)\n",
    "        timestep = lmp.style.timestep(units)\n",
    "        info = '\\n'.join([\n",
    "            'velocity move create %f %i mom yes rot yes dist gaussian' % (start_temp, randomseed),\n",
    "            'fix nvt all nvt temp %f %f %f' % (temperature, temperature,\n",
    "                                               tdamp),\n",
    "            'timestep %f' % (timestep),\n",
    "            'thermo %i' % (runsteps),\n",
    "            'run %i' % (runsteps),\n",
    "            ])\n",
    "    \n",
    "    return info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Run calculation function(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_dict = dislocationmonopole(lammps_command, ucell, potential, C,\n",
    "                                   burgers, ξ_uvw, slip_hkl,\n",
    "                                   mpi_command = mpi_command,\n",
    "                                   m = m,\n",
    "                                   n = n,\n",
    "                                   shift = shift,\n",
    "                                   shiftscale = shiftscale,\n",
    "                                   shiftindex = shiftindex,\n",
    "                                   sizemults = sizemults,\n",
    "                                   amin = amin,\n",
    "                                   bmin = bmin,\n",
    "                                   cmin = cmin,\n",
    "                                   etol = energytolerance,\n",
    "                                   ftol = forcetolerance,\n",
    "                                   maxiter = maxiterations,\n",
    "                                   maxeval = maxevaluations,\n",
    "                                   dmax = maxatommotion,\n",
    "                                   annealtemp = annealtemperature, \n",
    "                                   annealsteps = annealsteps,\n",
    "                                   randomseed = randomseed, \n",
    "                                   boundaryshape = boundaryshape,\n",
    "                                   boundarywidth = boundarywidth,\n",
    "                                   boundaryscale = boundaryscale)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['dumpfile_base', 'symbols_base', 'dislocation', 'dumpfile_disl', 'symbols_disl', 'E_total_disl'])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results_dict.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Report results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.1 Define units for outputting values\n",
    "\n",
    "- __length_unit__ is the unit of length to display results in.\n",
    "- __energy_unit__ is the unit of energy to display cohesive energies in.\n",
    "- __energy_per_length_unit__ is the energy per length to report the pre logarithmic energy factor in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "length_unit = 'angstrom'\n",
    "energy_unit = 'eV'\n",
    "pressure_unit = 'GPa'\n",
    "energy_per_length_unit = energy_unit+'/'+length_unit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.2 Print Stroh method parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pre-ln factor (alpha) = 0.38115202169620965 eV/angstrom\n",
      "K_tensor =\n",
      "[[ 81.07544348   0.         -12.76027037]\n",
      " [  0.         123.86918865   0.        ]\n",
      " [-12.76027037   0.         127.31323026]] GPa\n"
     ]
    }
   ],
   "source": [
    "print('pre-ln factor (alpha) =', uc.get_in_units(results_dict['dislocation'].dislsol.preln, energy_per_length_unit), energy_per_length_unit)\n",
    "print('K_tensor =')\n",
    "print(uc.get_in_units(results_dict['dislocation'].dislsol.K_tensor, pressure_unit), pressure_unit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.3 List dump files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Perfect system is saved as base.dump\n",
      "Defect system is saved as disl.dump\n"
     ]
    }
   ],
   "source": [
    "print('Perfect system is saved as', results_dict['dumpfile_base'])\n",
    "print('Defect system is saved as',  results_dict['dumpfile_disl'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
